
package com.jstarcraft.ai.jsat.math.integration;

import static java.lang.Math.abs;
import static java.lang.Math.pow;

import com.jstarcraft.ai.jsat.math.Function1D;

/**
 * This class provides an implementation of the Romberg method for numerically
 * computing an integral
 *
 * @author Edward Raff
 */
public class Romberg {

    /**
     * Numerically computes the integral of the given function
     *
     * @param f the function to integrate
     * @param a the lower limit of the integral
     * @param b the upper limit of the integral
     * @return an approximation of the integral of &int;<sub>a</sub><sup>b</sup>f(x)
     *         , dx
     */
    public static double romb(Function1D f, double a, double b) {
        return romb(f, a, b, 20);
    }

    /**
     * Numerically computes the integral of the given function
     *
     * @param f   the function to integrate
     * @param a   the lower limit of the integral
     * @param b   the upper limit of the integral
     * @param max the maximum number of extrapolation steps to perform.
     *            &int;<sub>a</sub><sup>b</sup>f(x) , dx
     */
    public static double romb(Function1D f, double a, double b, int max) {
        // see http://en.wikipedia.org/wiki/Romberg's_method

        max += 1;
        double[] s = new double[max];// first index will not be used
        double var = 0;// var is used to hold the value R(n-1,m-1), from the previous row so that 2
                       // arrays are not needed
        double lastVal = Double.NEGATIVE_INFINITY;

        for (int k = 1; k < max; k++) {
            for (int i = 1; i <= k; i++) {
                if (i == 1) {
                    var = s[i];
                    s[i] = Trapezoidal.trapz(f, a, b, (int) pow(2, k - 1));
                } else {
                    s[k] = (pow(4, i - 1) * s[i - 1] - var) / (pow(4, i - 1) - 1);
                    var = s[i];
                    s[i] = s[k];
                }
            }

            if (abs(lastVal - s[k]) < 1e-15)// there is only approximatly 15.955 accurate decimal digits in a double, this
                                            // is as close as we will get
                return s[k];
            else
                lastVal = s[k];
        }

        return s[max - 1];
    }
}
